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ABSTRACT 

Context. The direct detection of exoplanets with high-contrast imaging requires advanced data processing methods to disentangle 
potential planetary signals from bright quasi-static speckles. Among them, angular differential imaging (ADI) permits potential plan¬ 
etary signals with a known rotation rate to be separated from instrumental speckles that are either statics or slowly variable. The 
method presented in this paper, called ANDROMEDA for ANgular Differential OptiMal Exoplanet Detection Algorithm is based on 
a maximum likelihood approach to ADI and is used to estimate the position and the flux of any point source present in the field of 
view. 

Aims. In order to optimize and experimentally validate this previously proposed method, we applied ANDROMEDA to real 
VLT/NaCo data. In addition to its pure detection capability, we investigated the possibility of defining simple and efficient crite¬ 
ria for automatic point source extraction able to support the processing of large surveys. 

Methods. To assess the performance of the method, we applied ANDROMEDA on VLT/NaCo data of TYC-8979-1683-1 which is 
surrounded by numerous bright stars and on which we added synthetic planets of known position and flux in the field. In order to 
accommodate the real data properties, it was necessary to develop additional pre-processing and post-processing steps to the initially 
proposed algorithm. We then investigated its skill in the challenging case of a well-known target, f} Pictoris, whose companion is 
close to the detection limit and we compared our results to those obtained by another method based on principal component analysis 
(PGA). 

Results. Application on VLT/NaCo data demonstrates the ability of ANDROMEDA to automatically detect and characterize point 
sources present in the image field. We end up with a robust method bringing consistent results with a sensitivity similar to the recently 
published algorithms, with only two parameters to be fine tuned. Moreover, the companion flux estimates are not biased by the 
algorithm parameters and do not require a posteriori corrections. 

Conclusions. ANDROMEDA is an attractive alternative to current standard image processing methods that can be readily applied to 
on-sky data. 

Key words. Methods: data analysis - Techniques: high angular resolution - Techniques: image processing - Planetary systems - 
Planets and satellites: detection 


1. Introduction 

Among the thousands of exoplanets known today, less than 2% 
have been detected by direct imaging. The methods most of¬ 
ten used at the present date are indirect, such as Doppler spec¬ 
troscopy (Marcy & Butler 1993; Walker 1995), transit photome¬ 
try (Rosenblatt 1971; Borucki et al. 1985), or microlensing (Cas- 
san et al. 2012), that have permitted the population of inner plan¬ 
ets (semi-major axis < 5 AU) of sub-Jovian mass (< 2.10^*^ g) to 
be characterized (Udry & Santos 2007; Schneider et al. 2011). 
Direct imaging is motivated by its ability to detect massive 
young exoplanets on wider orbits and also to obtain spectra that 
provide information about their atmospheric composition. Thus, 
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direct imaging, complementary with other indirect techniques, 
increases our statistics and understanding of exoplanets. 

Very few exoplanets have been imaged up to the present day, 
and this is essentially due to the difficulty in achieving high con¬ 
trast at small angular separation, as most giant exoplanets lie 
within 1" from their host star with a flux ratio of about 10^ to 
10* (depending on mass and age). In order to obtain a high res¬ 
olution, observation with very large ground-based telescopes is 
needed, but the resolution is then limited by atmospheric tur¬ 
bulence. This effect can be overcome by the use of an adaptive 
optics (AO) system, which performs a real-time correction of 
the incoming distorted wavefront. The second step is to han¬ 
dle the very high contrast that exists between the planet and its 
host star. The use of a coronagraph greatly helps by removing 
a large part of the coherent light from the star while preserv¬ 
ing its close environment. However, the optics are not perfect, 
and in long exposures the remaining aberrations are responsi- 
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ble for speckles in the image that vary from less than minutes 
to hours (Hinkley et al. 2007); the speckles are too slow to be 
averaged into a smooth halo and too fast to be calibrated and re¬ 
moved a posteriori. These quasi-static speckles make exoplanet 
detection confusing since they are of the same angular size as the 
expected planetary signal and are often brighter. Consequently, 
at this stage it is necessary to apply elaborated image process¬ 
ing methods to disentangle the signal of an exoplanet from the 
remaining speckle noise. 

The image processing methods that are the most widely used 
by astronomers today are based on the angular differential imag¬ 
ing (ADI) concept (Marois et al. 2006). This technique relies on 
observations made in pupil tracking mode that fixes the pupil 
while the image rotates with time. Several algorithms are used 
to build a reference point spread function (PSF) that should re¬ 
produce as accurately as possible the speckle pattern to be sub¬ 
tracted, but not the rotating signature of a potential astronomi¬ 
cal source around the on-axis star. Among them. Classical ADI 
(cADI) uses a median image to build the reference PSF, and 
LOCI (Locally Optimized Combination of Images) (Lafreniere 
et al. 2007) uses a linear combination of images to build it and 
locally optimizes the speckle subtraction. The principal compo¬ 
nent analysis (PCA), in the form of PynPoint (Amara & Quanz 
2012) or KLIP (Soummer et al. 2012), uses an expansion of 
eigenimages to build the reference PSF. Numerous variations of 
these methods are studied and are used to improve the speckle 
subtraction. 

The ANgular DiffeRential OptiMal Exoplanet Detection Al¬ 
gorithm (ANDROMEDA) methodpresented in this paper is an 
advanced way of exploiting the ADI technique in the framework 
of inverse problems. ANDROMEDA uses a maximum likeli¬ 
hood estimation to detect planetary signals and retrieve the pro¬ 
jected distance between the planet and the star, as well as the 
contrast between the two components. Erom these estimations, 
and knowing the age of the host star, sophisticated models (e.g., 
Allard et al. (2012); Spiegel &. Burrows (2012); Marleau &. Cum- 
ming (2014) for young stars) provide estimations of the compan¬ 
ions’ mass and surface temperature. These results are useful to 
constrain the current models of star and planetary formation and 
their evolution. 

In this paper we first present and explain in detail the differ¬ 
ent steps performed by the algorithm ANDROMEDA to properly 
process real images and we discuss the hypothesis made to check 
whether it is consistent for real data (Sect. 2). We then explain 
how to use the output provided by ANDROMEDA to automat¬ 
ically detect the point sources in the images and to extract both 
their flux and position (Sect. 3). The next section is an applica¬ 
tion of the full procedure to real VLT/NaCo data, which con¬ 
sists of a bright star (TYC-8979-1683-1) suiTounded by numer¬ 
ous background stars acting as point sources to be detected, and 
in which we injected synthetic companions in order to assess the 
performance of ANDROMEDA (Sect. 4). The following section 
shows results using VLT/NaCo data of the well-known case of 
Beta-Pictoris, imaged under fair conditions and using an AGPM 
coronagraph (Mawet et al. 2005), thus approaching the quality of 
future high-resolution and high-contrast imaging systems (Absil 
et al. 2013). This last part (Sect. 5) compares the performance 
obtained by using ANDROMEDA with another commonly used 
method, the PCA. 

2. ANDROMEDA’S principle 

ANDROMEDA uses a statistical approach to discriminate plan¬ 
etary signals from the remaining speckles. The first step is to per¬ 


form ADI in order to create differential images in which speck¬ 
les that are stable enough to be distinguished from a rotating 
companion signature are removed. If a companion is present, a 
specific signature appears at its location that we can model ac¬ 
cording to its flux and initial position parameters. The second 
and original step of the method is to search for such a signature 
in the differential images by using a maximum likelihood esti¬ 
mation (MLE) of its position and intensity. This estimation is 
made from a consistent model of the data set, given a statistical 
model of its noise distribution as first presented in Mugnier et al. 
(2009), it is similar to the approach proposed for the processing 
of DARWIN/TPE-I data in Thiebaut & Mugnier (2006) and to 
the works by Smith et al. (2009) for a perfectly fixed speckle 
field. 

The purpose of this paper is to apply ANDROMEDA to real 
data from the VLT/NaCo instrument (Lenzen et al. 2003; Rous- 
set et al. 2003). This implied adapting the method described in 
Mugnier et al. (2009). The different processing steps performed 
by the upgraded version of ANDROMEDA presented and used 
in this paper are the following. Section 2.1 describes how and 
why the artifacts with low spatial frequency are removed from 
the raw images in a pre-processing step. Section 2.2 describes 
the angular image difference imaging that is then performed in 
order to suppress the vast majority of the starlight from the im¬ 
ages. Section 2.3 defines the mathematical model corresponding 
to the resulting differential images along with their dependency 
with the sought positions and amplitudes of the potential com¬ 
panions. Section 2.4 explains how the problem is inverted be¬ 
tween the differential image and its model by using a maximum 
likelihood which analytically finds an estimation of both the po¬ 
sition and flux of the potential companions. Einally, Sect. 2.5 
describes the simple and effective post-processing that is per¬ 
formed on the results in order to compensate for the deviation 
from the noise model (assumed to be white and Gaussian in the 
differential images) that appears in practice on real data. 


2.1. Pre-processing: Filtering out low-frequency artifacts 

In both the raw and the reduced images there are some strong 
inhomogeneities of low spatial frequency that are not stable and 
thus disturb planetary detection. These large-scale structures in 
the images, which vary slowly from one frame to another, might 
originate from temporal variation of the residual turbulence. As 
they cannot be fully subtracted via the ADI process and they are 
not included in the image model used to perform the MLE, these 
low frequencies must be removed. 

A pre-processing of the data has thus been introduced in AN¬ 
DROMEDA to make the detection of point sources easier within 
the stellar halo. This pre-processing consists in eliminating these 
disturbing structures by applying a high-pass filter in the Eourier 
space. Half of a Hann function has been chosen to build this fil¬ 
ter’s profile in order to avoid Gibbs effects, due to discontinuity, 
that appear when going from the Eourier to the real space. This 
filter attenuates the spatial frequencies lower than a chosen cut¬ 
off frequency defined as fc - F ■ fy, where is the Nyquist 
frequency and F is a factor in the range 0-1 (hereafter called 
the filtering fraction). The parameter F is user-defined and must 
be chosen carefully to efficiently remove the low-frequency ar¬ 
tifacts while preserving most of the energy of the point source 
signals in the images. In order to quantify the signal loss due to 
this filtering, we took an Airy pattern and applied the filtering to 
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it. The energy loss is then quantified as follows: 

_ j (airy filtered) 

^loss “ip, 2 ' 

J (airyjton-filtered) 

The obtained results are gathered in a graph showing the energy 
loss as a function of the filtering fraction (Fig. 1). 
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Fig. 1. Energy loss as a function of the filtering fraction when filtering a 
simple Airy pattern (sampled at the Nyquist rate) with a high-pass Han¬ 
ning Fourier filter (solid line). Energy loss for a usual filtering fraction 
of 1/4 is shown with the red dashed line. 

Comparisons and tests on simulation and real data show that 
even though this procedure slightly reduces the intensity of the 
planetary signal on individual frames, the filtering is worthwhile 
for planet detection and does not introduce any etTor on the esti¬ 
mated flux. In practice, a filtering fraction of about one fourth is 
recommended when dealing with VLT/NaCo data, leading to an 
energy loss of about 18% (see Fig. 1). 

2.2. Construction of differential images 

2.2.1. Angular differential imaging (ADI) 

The ADI technique consists in taking advantage of the alt- 
azimuthal mount of the telescope in use: the pupil and the field 
rotate at different deterministic velocities during the observation. 
In order to perform ADI, one can choose to fix the pupil (by us¬ 
ing a pupil derotator if needed) to stabilize the aberrations in the 
image, while the observed field rotates as the parallactic angle. 
Therefore, these images can be later subtracted from each other 
in an optimized way to reduce as much as possible the speckle 
noise in the resulting differential image. This method was first 
introduced by Marois et al. (2006) and is the basis of most im¬ 
age processing methods dedicated to exoplanet detection today. 

In order to choose the two images to be subtracted, a com¬ 
promise has to be made so that (i) the turbulence strength and the 
quasi-static abetTations do not have time to evolve significantly 
between the two images, and (ii) the object of interest has rotated 
enough in the field between the two images in order not to can¬ 
cel its own signal during the subtraction. Assuming that the most 
correlated data are the closest in time, we chose two images that 
are as close as possible in time but within a limitation to avoid 
the self-subtraction of the companion. This limitation is the min¬ 
imum distance between the position of the planetary signal in the 
two images (dmin), which is a key user-defined parameter whose 


value must be adjusted according to the speckle lifetime and in¬ 
strumental stability (typically l.QA/D, see Table 2). As the dis¬ 
placement between two images is dependent upon the distance 
between the star and the object of interest, the ADI is performed 
over concentric annuli of constant width dr to accommodate this 
dependency, which is another user-defined parameter (typically 
l.QA/D, see Table 2). For each of these annuli, each image of the 
cube is subtracted from the first following image (or previous 
image if needed) that meets condition (ii). If a couple is found, 
the subtraction results in the creation of a differential image, in¬ 
dexed k, where k e [1; A*;], with Nk the total number of couples 
found within the cube at the regarded distance from the star. If 
for a given frame, no other frame satisfies condition (ii), no pair 
can be built with this frame which is therefore not used. Thus Nk 
is less than or equal to the total number of images in the cube. 

As the two images to be subtracted are likely not to have the 
same intensity distribution, one of the image needs to be adjusted 
with respect to the other via a scaling factor, noted jk- These in¬ 
tensity variations are primarily caused by variations of the seeing 
in the course of the observation. Each differential image A(r, k) 
resulting from the subtraction of a couple k of images i(r, t) is 
then of the form 

A(r, k) = ik(r, fi) - ^ 4(r, tz) . (2) 

It can be seen in both the simulations and the real data that the 
flux difference depends on the distance to the star. As a conse¬ 
quence, the scaling factor is computed for each couple of annuli, 
over an optimization area that is slightly wider than the effec¬ 
tive subtraction area in order to avoid discontinuities between 
adjacent annuli (Cornia et al. 2010). The ratio optimization to 
subtraction area (R^) is a third user-defined parameter, set con¬ 
stant for all annuli (typically 2.0, see Table 2). 

Several considerations guided us in the design of the best 
computation of this scaling factor. Because the high-pass filtered 
images have a zero mean, it is not possible to compute the scaling 
factor jk with a simple ratio of the total intensities in each image. 
Instead the scaling factor is better estimated by a least-squares 
fit that minimizes ||4(fi)(r) - 7k 4(/2)(r)lP (Cornia et al. 2010). 
In the following, we refer to this method as the LS optimization. 
We note that optimization depending on position in the field has 
been implemented in other methods of exoplanetary detection, 
as in Marois et al. (2006); Lafreniere et al. (2007). Another way 
to subtract the images is also presented here, after taking two 
aspects into consideration: first the PSF may not vary linearly 
and if there is a sudden evolution in the turbulence strength, the 
image can be better fitted by an affine law. Thus, the differential 
image can be constructed by 

A(r, k) = 4(,2)(r) - jk 4pi)(r) - y'k- (3) 

Second, the least-squares is sensitive to outliers and if, from one 
image to the other, a speckle or a planetary signal intensity has 
varied significantly, their high signal is taken into account to 
compute the scaling factor, which results in a flux offset in the 
differential image. In order to avoid the latter bias, we addition¬ 
ally implemented a robust fit, using a Li norm that can be chosen 
instead of the L 2 norm. The scaling factors jk and are then 
computed by minimizing: ||4(r) - jk ■ 4(r) - T^IIli • In case there 
are such outliers corresponding to high flux and variable signals 
(either speckles or planetary signals), this way of combining the 
image pair - hereafter called LI-affine optimization - has proven 
equally or more efficient at reducing the speckle noise on both 
simulated and real data, while preserving the planetary signal, 
if any. This LI-Affine optimization is particularly good at very 
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close separation where the intensity of the numerous speckles 
varies a lot from one image of the couple to the other due to the 
poor rotation field velocity. 

After this operation, the differential images A(r, k) are ob¬ 
tained for every distance to the star. We note that even if the 
image couples are determined for specific annuli and the opti¬ 
mizations are computed over annular areas including these an¬ 
nuli, the whole images are actually subtracted one from another 
to build the differential images. By minimizing the residuals in 
these data, the noise has been partly whitened spatially (Cor- 
nia et al. 2010; Mawet et al. 2014). Since ANDROMEDA relies 
on solving an inverse problem, it is also useful to estimate the 
variance map cr^(r) of the differential images. Currently it is es¬ 
timated empirically over the time dimension k as a function of 
the position r; 

crlir) = <[A(r, k)]\ - [<A(r, k)),]\ (4) 

We note that when the observer has other independent ways to 
estimate the variance map, in particular taking the temporal vari¬ 
ations into account, it could be introduced here. In the following, 
it is nevertheless assumed that the noise variance is stationary in 
time and the definition of the variance map given above is used, 
which is satisfactory for our study. 

2.2.2. Other possible combinations of images 

Because ANDROMEDA was developed with the objective of 
being used to process large amount of survey data from second- 
generation instruments, it was proposed with a simple ADI sub¬ 
traction to whiten the noise. This allows the algorithm to run 
quickly (in particular, it is possible to parallelize the computa¬ 
tion of each operation), and to have a very low number of user- 
parameters to be adjusted. However, many other approaches can 
be considered to build up differential images that could also be 
included in the ANDROMEDA scheme, as long as the impact 
on the potential companion signature can be evaluated and in¬ 
serted in the model. Eor instance, it would be possible to merge 
this method with other subtraction algorithms, such as LOCI or 
PCA, that are known to be able to reduce the speckle noise very 
efficiently. 

Other techniques can be used to whiten the noise: instead 
of exploiting only the temporal information, ANDROMEDA 
could deal with spectral or polarization redundancy, for instance, 
and perform, respectively, spectral differential imaging (SDI) 
(Racine et al. 1999; Marois et al. 2000) or polarimetric differ¬ 
ential imaging (PDI) (Kuhn et al. 2001) in addition to the ADI. 
An optional SDI for ANDROMEDA has been implemented and 
discussed in Cornia et al. (2010) and simple tests have been led 
on images taken either with a dual band imager (DBI) or an 
integral field spectrograph (lES), such as IRDIS or lES on the 
VLT/SPHERE instrument, but this work is beyond the scope of 
this paper. 

2.3. Building a model of the differential images 

Once the differential images have been created, evidence of the 
presence of a planetary signal in the original image is sought in¬ 
side all of these differential images. After performing the ADI, 
if a companion is indeed present in the field of view, a pecu¬ 
liar signature appears in the differential images. This so-called 
planet signature is naturally the difference of two planetary sig¬ 
nals separated by at least dmin. Eigure 2 illustrates the shape of 
the expected planet signature obtained either with or without 



Fig. 2. Illustration of the planet signatures obtained when performing 
the ADI, setting d„„>, = 0.5A/D. These have been obtained by simulating 
two identical noiseless planetary PSF, separated by 6min and then sub¬ 
tracted one from the other. Left: Planet signature obtained without high- 
pass filtering of the raw data. Right: Planet signature obtained when the 
raw data have been filtered of their low frequencies {F = 1/4). This time 
secondary opposite lobes surround the main lobe. 

high-pass filtering of the raw data. In order to properly estimate 
the position and flux of a potential exoplanet using a maximum 
likelihood approach, we now define a model for these differen¬ 
tial images. Assuming that the stellar halo is fully suppressed by 
the subtractions and that a planet is present in the field of view, 
each differential image can be modeled as 

A(r, k) - a ■ p(r, k; Tq) + n(r, k) , (5) 

where the planet’s flux a and the position of the planet on the first 
image of the sequence (initial position) Tq are the two unknown 
parameters, n(r, k) denotes the residual noise, and p(r, k; Tq) is 
the planet signature. 

The model can be built by placing two PSEs, one positive 
and one negative, properly positioned in the field of view to cor¬ 
respond to the field rotation at times fi and t 2 of a certain couple 
k. We assume here that the exposure time is short enough that 
the companion does not undergo any azimuthal smearing, but it 
could be added if the smearing becomes significant. As we have 
only one representative PSE for all the images of the cube, we 
also assume that the PSE is known and does not vary with time. 
As was true for other ADI approaches, if the instrument PSE 
does vary during the sequence, it directly affects the companion 
flux estimate. The main limitation here is related to monitoring 
the PSE Strehl ratio along the sequence. If it were known, this 
knowledge could be included to rescale the searched companion 
signature in ANDROMEDA within the succesive differential im¬ 
ages. In the current paper, we also work under the hypothesis that 
the companion PSE core is temporally stable. In this case, the 
scaling factor used to build up the differential images, which 
is related to the variability of the star halo, must be taken into 
account to build the model of the planet signature, 

p(r, k; ro) = h(r - r'(fi, ro)) - yk h(r - r'(f 2 , ro)), (6) 

where h is the PSE of the system that can be estimated by provid¬ 
ing an empirical PSE to the algorithm simply obtained by imag¬ 
ing the star (unsaturated image or without coronagraph). This 
PSE is shifted at a position r'(fj, ro) that the planet will have at a 
time tx if its position in the first frame is Tq. If, however, it turns 
out that the intensity of the planetary signal and the PSE vary the 
same way in every image, the scaling factor should not be taken 
into account to build this model. 

If another kind of method is used to perform the speckle sub¬ 
traction, as suggested in Sect. 2.2.2, the model must be modified 
according to the subtraction performed. Eor instance, if the pre¬ 
viously mentioned methods LOCI or PCA are used, the coef¬ 
ficients found to perform the subtraction must be recorded for 
each zone in order to build a consistent model. 
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Since ANDROMEDA models differential images through 
Eq. (5), the filtering procedure must also be applied to the ref¬ 
erence PSE used to build the planet signature. In this way, the 
model remains consistent with the data and ANDROMEDA 
tracks the proper filtered planet signature in the filtered data se¬ 
quence. 


2.4. Estimation of position and fiux via maximum iikeiihood 

ANDROMEDA’S goal is to estimate the position and the flux of 
a potential planet orbiting around the target star using a proba¬ 
bilistic approach. This estimation is based on a MLE under the 
hypothesis that the residual noise inside the differential images 
is a Gaussian and white. Ear from the star, this hypothesis is ful¬ 
filled if the star PSE does not vary during the exposure because 
there is only photon and detector noise in the images (Mugnier 
et al. 2004). Closer to the star, quasi-static speckles are predom¬ 
inant and they follow modified Rician statistics (MR) (Good¬ 
man 1985; Soummer & Aime 2004; Eitzgerald & Graham 2006), 
which are neither Gaussian nor totally white. 

However, ANDROMEDA deals with high-pass filtered im¬ 
age differences. Eirst, it appears that the filtering process of the 
raw images increases the Gaussianity of the residual distribu¬ 
tion in the images since it removes the slowly evolving struc¬ 
tures while preserving the high spatial frequencies that evolve 
quickly. As the remaining noise inside an annulus is with a good 
approximation independant, according to the central limit theo¬ 
rem, its distribution should follow a normal distribution. Second, 
the speckle noise is temporally correlated from one frame to an¬ 
other. By performing the difference between close images, only 
the variation between the speckle fields in the two frames will 
remain in the differential image. That is to say the created dif¬ 
ferential image will contain the time derivative of the speckle 
and the power spectral density (PSD) of a derivative tends to¬ 
ward a constant so the differential imaging process “whitens” the 
noise. This whitening is a well-known consequence of the ADI 
(Marois et al. 2008), which is true mostly because we consider 
the reference frame of the sky. This property is verified empiri¬ 
cally on real NaCo data in Sect. 4.2.2. To conclude, the noise in 
the differential images created to perform the MLE is, to a good 
approximation, Gaussian and white, thanks to the filtering and 
differential imaging performed. At the end of this section there 
is further discussion of this hypothesis. 

Under the assumption that the residual noise in the differ¬ 
ential images is white (in time and space), Gaussian, and non¬ 
stationary, the likelihood can be written 


L(r„,a) ocexp{-i^^ 


|A(r, k)-a p(r, k; ro)F 
mlir) 


(7) 


Eor each annular zone, at a certain distance from the star, this 
likelihood is computed for any possible initial position Tq: the 
sum is made over all the couples found for this distance from 
the star and over all the pixels in the image field. 

ANDROMEDA analyzes the sequence of differential images 
A(r, k) by seeking the optimal values (fo, a) that maximize the 
logarithm of the likelihood. 


7(ro, a) A In L(ro, a) oc - ^ 

r.k 


|A(r,k) - a ■ p(r, k; ro)r 
2o-2(r) 


( 8 ) 


the residual noise. We note that the definition of the variance cr^ 
at Eq. 4 means that the weight is lower when there is a planet 
in the studied differential image (and the noise is more overesti¬ 
mated to a higher degree when the planet is brighter and at closer 
separation). This bias induces an overestimation of the error, but 
this does not affect the estimation itself. 

The optimal flux value a for each possible initial position of 
the planet is easily computable analytically and is given by 


'Lr.k pjr, k; ro)A(r, k)lo\{r) 
ro)/cr2 (r) 


(9) 


This equation shows that the estimated flux can be regarded as 
a cross-correlation between the planet signature and the differ¬ 
ential image, weighted by the noise variance averaged on every 
differential image, the denominator being a normalization con¬ 
stant. In the ANDROMEDA code, after a processing per annuli, 
the data are combined in a single 2D map called a flux map, giv¬ 
ing at each pixel the flux of the object, if the object has this pixel 
position. 

By inserting this expression of a(ro) into the metric J, it is 
possible to obtain an expression for the log-likelihood that de¬ 
pends from now on only on the initial planet position Tq, 


7'(ro) A 7[ro, a(ro)] 


(2r,i k\ ro)A(r, k)lcr\{r)f 

'Lr.k E^(r, k ro)/o-2(r) 


+ C, (10) 


where C is a constant. The criterion 7'(ro) is a measure of the 
probability that there is a point source at position ro on the first 
image of the sequence. This formula is easily interpretable by 
saying that the planet has a high probability of being found at the 
position where the correlation between the model and the differ¬ 
ential image is the closest to one. In practice /'(ro) is computed 
for each possible initial position of the planet on a grid chosen 
as the original pixel grid of the images. The results are then a 
so-called likelihood map which has higher values at positions 
where the presence of a companion is more probable. 

Once the likelihood and flux maps are computed, one pur¬ 
pose is to link the intensity of a likelihood peak to a probabil¬ 
ity of false alarm in order to assess whether it is indeed a true 
planetary signal. One way to proceed is to compute the standard 
deviation of the estimated flux, cr[a(ro)], for each possible initial 
position, which describes how the noise propagates from the dif¬ 
ferential image toward the flux map. The standard deviation of 
the estimated flux is given by 


o-2[a(ro)] = 


^/(r, k; ro)/cr^(r) 

r.k 


1 


( 11 ) 


It is then possible to define the signal-to-noise ratio of the esti¬ 
mated flux, S/N, as 


S/N(ro) A 


Q(ro) 

cr[fl(ro)] ■ 


( 12 ) 


The so-called detection limit, that is to say the faintest detectable 
companion flux as a function of the separation from the star, at 
chosen threshold Ncr is thus given by A x cr[a(ro)]. In other 
words, the standard deviation of the flux map is the detection 
limit at Icr. 

It is straightforward to show that the S/N can be rewritten as 


Equation (8) shows that maximizing the log-likelihood is equiv¬ 
alent to minimizing the sum of squared residuals between the 
differential image and the model, weighted by the variance of 


S/N(ro) = 


'Lr.k E(r, k-, ro)A(r, k)/cr2 (r) 
^I'Lr.k pHr, k; ro)/cr2(r) 


(13) 
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Showing that the log-likelihood is simply the square of the S/N. 
This means that maximizing the likelihood map is equivalent to 
maximizing the S/N map to find the most-likely position of a po¬ 
tential companion. The S/N thus gives a physical meaning to the 
log-likelihood J'{ro). Moreover, thresholding the S/N to zero is 
equivalent to computing a(ro) under the constraint that it is non¬ 
negative, and the log-likelihood then obtained through Eq 13 in¬ 
corporates this positivity constraint (see Mugnier et al. 2009). In 
summary, the likelihood and S/N maps contain the same infor¬ 
mation, and we chose to use the S/N map to perform detection 
since it directly yields the statistical significance of each detec¬ 
tion. 

It is then possible to define a threshold r, to be applied to 
the S/N map, which corresponds to a certain confidence level of 
the detection. Assuming that there is no planet anywhere in the 
images (hypothesis Hq) our model shows that the differential im¬ 
ages equal the residual noise. If our hypothesis were strictly full- 
filled, the residuals in the S/N map would have the same statistics 
as the noise in the differential images (because each pixel of the 
S/N map is a linear combination of pixels from the differential 
images (via Eq. (9) and (12)), which is assumed to be white and 
Gaussian, with a probability density function (PDE), ps/^ that 
is normal, zero-mean, and of unit variance. Under the assump¬ 
tions on the noise, the probability of false alarm (PEA), defined 
as the probability that the S/N shows a signal above the chosen 
threshold r under the hypothesis Hq, writes 

PEA(t) = PriS/N > t)|h„ 

1 r’’ 

= 1-— I exp(-T'^/2) dr' 

y7jt J-oc 

^l-erf(T), (14) 

where erf is the so-called error function. The applicability and 
relevance of such a threshold t is important as further discussed 
in Sect. 2.5 and when applying ANDROMEDA over real data in 
Sect. 4. 

However, these considerations are for the ideal case of a 
white, Gaussian, and non-stationary noise in the differential im¬ 
ages, for which the connection between PEA and threshold is 
perfect. Eor real data, even if the filtering permits an increase in 
the Gaussianity of the noise and allows the differential imaging 
to remove an important fraction of its correlated component, the 
actual noise distribution still slightly deviates from our assump¬ 
tions. The real noise distribution is closer to a MR distribution, 
which increases the PEA. Because of its structure, the MR distri¬ 
bution evolves radially; the quasi-static speckles have the same 
mean and variance only for a specific annulus. Close to the star, 
the speckle noise is dominant (and the noise is more correlated), 
whereas far from the star the photon noise, the readout noise, and 
the dark current prevail. The threshold should thus be increased 
according to the separation from the star where the companion 
is sought (Marois et al. 2008). 

Another effect that is worth pointing out here is that for an¬ 
nuli close to the star, the field has not rotated as much as at large 
separations. At small separations, not only are there a small num¬ 
ber of pixels contained in the annulus, there are also few images 
to be subtracted. These two effects reduce the number of degrees 
of freedom (in ANDROMEDA, the number of points inside the 
annulus, r, plus the number of ADI couples, k) and corrupt the 
statistics; we are in a small sample statistics regime. In that case, 
Mawet et al. (2014) have shown, using a frequentist approach, 
that the PEA is underestimated to an even greater extent. To cor¬ 
rect for this bias, close to the star, the noise statistics no longer 
follow a MR distribution, but a robust Student’s t-distribution. 


which is valid only if the sample is independent and identically 
distributed (which is the case if processes to whiten the noise, 
such as differential imaging, have been applied to the images). 
In practice, when correcting for this bias, we calculate the equiv¬ 
alent threshold that should be applied to provide the same PEA 
that Gaussian statistics would give at the actual chosen thresh¬ 
old. 

As a consequence, the PDE chosen in Eq. (14) is no longer 
the normal law, N, and should be replaced by the appropriate 
law; a MR PDE at large separation (Marois et al. 2008) and a Stu¬ 
dent t-distribution PDE at small separation (Mawet et al. 2014). 
The detection limit calculation should then be modulated to take 
into account these two considerations. The next section is about 
how to take the deviation into account from our model, in order 
to set a constant threshold throughout the whole S/N map. Set¬ 
ting a constant threshold enables the probable point sources in 
the field to be automatically and systematically detected with a 
relevant confidence level. 

We note that the discrepancy between the model of noise 
made with respect to the real noise distribution does not affect 
the flux estimation in itself, but will affect the error on the flux 
estimation. 

2.5. Post-processing: normalization of the S/N map 

To complete the presentation of the method within this section, 
we can already indicate that the first tests performed on real data 
showed that the deviation from the noise model prevents one 
from applying a constant threshold to the S/N map and building 
an automatic procedure to detect companions. We can vizual- 
ize this deviation effect in Eig. 9-left, which notably shows a ra¬ 
dial dependancy so the threshold should indeed be modified as a 
function of the separation. 

This is not unexpected since if the noise model (white, Gaus¬ 
sian, and non-stationary) were consistent with the real values, the 
output S/N map from ANDROMEDA would have a zero mean 
and a standard deviation of one when, of course, excluding the 
pixels containing the signal from a companion. In this case, the 
threshold could be uniform all over the field. However, as ex¬ 
plained in the previous section, in real images, the noise clearly 
deviates from this model and the standard deviation of the resid¬ 
ual noise in the S/N map generally has a standard deviation larger 
than one which increases from the center to the edge and a mean 
of zero very far from the star that increases when going closer to 
the star. 

Knowing that we would like to obtain a S/N map that has a 
zero mean and a standard deviation of one, we propose a simple 
solution to preserving the PEA, which consists in normalizing 
radially the S/N map by its own empirical radial standard devi¬ 
ation. This operation removes the radial dependency and allows 
one to apply a constant threshold to the S/N map, regardless of 
the distance from the star. To normalize the S/N map, it is nec¬ 
essary to build an empirical radial profile of the standard devi¬ 
ation from the S/N map, without taking into account any peaks 
due to planetary signals. A method from Hoaglin et al. (1983) 
and Beers et al. (1990) to estimate this standard deviation, here¬ 
after called robust standard deviation, has been implemented in 
ANDROMEDA. This method consists in calculating the median 
absolute deviation divided by a normalization factor enabling a 
robust estimate of the standard deviation. The factor is chosen 
so that the regular standard deviation and the robust value are 
identical in the case of a Gaussian distribution. Because only the 
global trend is needed, the mean radial robust standard devia¬ 
tion profile is smoothed over a certain number of pixels, Nsmooth. 
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which is a user-defined parameter that has to be set empirically 
according to the pixel scale and the size of the images. Once 
a satisfactory profile that does not take into account the small 
variations from one annulus to its neighbor is obtained, the nor¬ 
malized S/N map is obtained by dividing the S/N map by this 
radial trend. 

This method is reproducible and independent of the data set. 
The post-processed S/N map thus gives the confidence level for 
a point source detection on each pixel of the image grid. Thresh¬ 
olding this map by a constant value throughout the field of view 
provides a list of positions where a probable companion can¬ 
didate is detected. As a result, in the following we work exclu¬ 
sively with the normalized S/N map and from now on “S/N map” 
stands for the normalized S/N map. Through Eq. 12, as the S/N 
map has been normalized, the standard deviation flux map must 
also be normalized in a coherent manner to obtain a consistent 
relation between the normalized S/N and the standard deviation 
of the flux. The normalization process does not significantly af¬ 
fect the position and flux estimations (the impact is lower than 
the given errors on these estimations) but normalizing too much 
(smoothing the profile too much) artificially lowers the error on 
the flux estimation. 

ANDROMEDA provides four 2D maps for analysis: the like¬ 
lihood map, the S/N map, the flux map, and the flux standard 
deviation map. The two most useful outputs for detection and 
characterization are the S/N map, in which we can look for plan¬ 
etary signals and estimate their most likely position ro, and the 
flux map, in which we can read the corresponding estimated flux 
a{ro). The map of the flux standard deviation is the detection 
limit at Icr. The following section deals with the analysis of the 
ANDROMEDA output to accurately compute the position and 
flux of potential companions present in the image field based on 
applications on real VLT/NaCo images to show concrete illus¬ 
trations. 

3. Extracting companion information from the 
ANDROMEDA outputs 

In this section we explain how to systematically and efficiently 
extract a list of likely companions with their subpixel position 
and estimated flux in terms of contrast. Once the S/N and flux 
maps are at hand, a complementary module is developed that 
automatically detects and analyzes the point sources present in 
the images. This module simply finds the signals above threshold 
in the S/N map and, knowing the expected shape of a planetary 
signal in the S/N map (Sect. 3.1), assesses their subpixel posi¬ 
tions and then estimates their corresponding flux thanks to the 
flux map (Sect. 3.2). Also, during the analyses, some tests are 
performed to discriminate probable planetary signals from arti¬ 
facts (Sect. 3.3). 

3.1. Resulting pattern from a planetary signal 

If there is indeed a planetary signal, a specific pattern appears 
in both the S/N and flux maps. This pattern is well known and 
derived from Eq. (10) and Eq. (13): under the hypothesis that 
there is no noise in the data, it can be seen essentially as the 
autocorrelation of the planet signature p{r, k; To) shown in Eig. 2, 
multiplied by the intensity of the candidate companion. Such a 
pattern is illustrated in Eig. 3, either with or without pre-filtering 
of the data. 

As expected, this pattern is made up of a main positive lobe 
surrounded by two negative lobes positioned perpendicular to 


the star-companion direction. The main lobe of the pattern is 
oval with its longest side along the radial direction with a typ¬ 
ical length of A/D, whereas the length of its perpendicular di¬ 
rection is dependent upon the chosen distance dmin- In the case 
without pre-filtering, positively thresholding the S/N map retains 
only the main lobe. On the other hand, the filtering induces the 
appearance of two tertiary sickle-shaped positive lobes surround¬ 
ing the main lobe as well as two tiny fourth negative lobes. Con¬ 
sequently, when applying a positive threshold, the tertiary lobes 
may remain in the thresholded map if the S/N is high enough. 



Fig. 3. Illustration of the pattern (autocorrelation of the planet signa¬ 
ture, Fig. 2) appearing in the maps as evidence of a planetary signal. 
Left: Pattern obtained without filtering the raw images. Right: Pattern 
obtained with high-pass filtering of the raw images (F = 1 /4). 


3.2. Identification, position, and flux retrieval of candidate 
companions 

The S/N map permits the identification of the most likely com¬ 
panions and the extraction of their subpixel position in the first 
image. The signals above the user-defined detection threshold 
are identified and classified according to decreasing peak S/N 
values. For each of them, a subwindow is extracted that fully 
encloses the central part of the candidate signal, as shown for 
instance in Fig. 4, and whose size is about 4A/D. To evaluate the 
subpixel position Fq of a candidate companion in a simple way, 
the planetary signal is interpolated to find the subpixel position 
of the maximum S/N, coinciding with the position of the plane¬ 
tary signal. The main lobe of the pattern expected from a planet 
is approximated by a Gaussian (this approximation is quite reli¬ 
able in the very central part and less valid at the edges), and the 
subpixel position is identified at the maximum position of the 2D 
Gaussian fit of this main lobe. The error on the position estima¬ 
tion is simply given by the Gaussian fit uncertainty taken at 3cr. 
If needed, this error is corrected for the surrounding noise when 
it is too intense and might bias the error estimation. 

To read the corresponding estimated flux a(ro), the same ap¬ 
proach is used by performing a 2D Gaussian fit of the plane¬ 
tary signal, but this time using the flux map provided by AN¬ 
DROMEDA. The extracted flux of the companion is then the 
value of the fitted Gaussian at the subpixel position, Fq, previ¬ 
ously retrieved on the S/N map. At this stage, ANDROMEDA 
returns the flux of the detection with respect to the flux of the 
input PSF. From this value, it is easy to derive the contrast be¬ 
tween the two components by using all the information about 
how the input PSF has been imaged, such as the integration time 
of the PSF with respect to the integration time of each image of 
the cube, any normalization factor, or the neutral density trans¬ 
mission if one was used to image the star. We emphasize here 
that this approach intrinsically takes all the pre-processing ap¬ 
plied to the data into account: the companion signature that is 
sought takes into account the pre-filtering, and then the subtrac¬ 
tion applied on image pairs separated by a small rotation differ¬ 
ence. As a consequence the estimated flux does not depend on 
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such parameters, and contrary to other approaches such as PCA 
or LOCI, there is no companion flux loss. The error at Icr on 
the flux estimation is directly read on the flux standard deviation 
map at the location ro retrieved earlier. The error bars provided 
by the algorithm are the corresponding values at 3cr. 

3.3. Rejection of known sources of artifacts 

The method described so far, based on thresholding the S/N map, 
does not prevent the detection of artifacts whose signals are nev¬ 
ertheless above the threshold. Any signal temporally varying in 
the images is not fully subtracted and may leave some high-level 
residuals in the differential images, which could mimic the ex¬ 
pected planet signature. This kind of artifact appears in partic¬ 
ular if the pupil held is not well stabilized or if the star is not 
well centered in the image during the observation. Therefore, it 
is needed to discriminate a posteriori probable planetary signals 
from remaining artifact signals. The rejection criterion setup is 
based on the morphological properties of the pattern expected 
from a planetary signal (Sect. 3.1). 

Three main sources of errors have been identified that can 
be used to efficiently reject these false detections. To illustrate 
these sources, practical examples are pointed out in the images of 
Fig. 4 that were produced thanks to data provided by NaCo (see 
Sect. 4 for more details on the data processed). The subimages 
cut in the S/N map and classified by decreasing S/N signals are 
displayed in Fig. 4-left. Below are listed the three sources of false 
detection that have been diagnosed: 

1. The expected pattern contains tertiary positive lobes induced 
by the filtering procedure (Fig. 3-right). If the S/N of the 
signal is high enough, it these tertiary lobes may be above the 
threshold too and may thus be regarded as detections. These 
detected tertiary lobes are indeed surrounding only the very 
high S/N signals, such as #1, #2, and #3 visible in Fig. 4- 
right via the dark blue symbol @. Three criteria are used 
to spot and reject these artifacts: their proximity to a high 
S/N signal, the S/N ratio between the regarded signal and its 
neighbor, and their peculiar sickle shape which prevents the 
2D Gaussian fit from converging. Such signals can be seen 
in Fig. 4-left such as #18 (surrounding signal #1) and #25 
(surrounding signal #3 - we note that some pixels are missing 
in the subimage because they had already been erased after 
detecting the contingent part of this tertiary lobe). 

2. The filtering may reveal the spider diffraction pattern in the 
S/N map that are not completely subtracted during the ADI 
(probably owing to the bad centering of the star or pupil 
tracking). Consequently, some signals can be found above 
the threshold that are actually artifacts due to this diffrac¬ 
tion pattern. However, if this is the case, the signal is quite 
elongated along the radial direction and it is thus possible to 
constrain the fit parameters to reject such signals. This kind 
of signal can be seen in Fig. 4-left, such as signal #31, which 
is clearly located in the spider diffraction pattern visible in 
Fig. 4-right. 

3. In areas where the speckle noise is quite high (close to the 
star), we might detect signals that originate from the speckle 
noise. This kind of signal is usually quite irregular compared 
to the expected pattern and can thus be rejected by constrain¬ 
ing the 2D Gaussian fit to be neither too wide nor too shifted 
in the subimage. 

Some other sources of instrumental error (hot pixels, etc.) 
can induce the detection of an actual artifact. To find all these ar¬ 
tifacts and reject them, three of the fit parameters are constrained 


(by using the IDL function mpfit2Dpeak.pro developped by 
Craig B. Markwardt, for instance): (i) The distance between the 
fit maximum and the center of the subwindow must be of less 
than a pixel, (ii) the FWHM of the fit must be consistent with the 
one expected from the pattern, and (iii) the long-axis orientation 
must be close to the radial direction. In the following, a detection 
whose fit has not converged or that does not respect one of the 
constraints listed above is referred as an “ill-fitted” companion 
and is therefore rejected. In Fig. 4, when the 2D Gaussian fit is 
unsuccessful one asterisk is placed above the subimage (such as 
the detections indexed #17 and #25). When the S/N decreases, 
the shape of the signals is less and less regular and is conse¬ 
quently often not well fitted. Indeed, the lower the S/N, the less 
probable it is for the signal to be real, hence the frequency of 
strangely shaped signals increases (e.g., we can see this effect 
in signal #31, which has a S/N value of 5.4cr and lower). If the 
signal found is too close to the edge of the image, so its signal 
might be incomplete, it is impossible to obtain a correct 2D fit 
and in that case three asterisks are marked above the subimage 
(e.g., signal #8). 

4. Application to NACO data: Test case with 
TYC-8979-1683-1 data 

In order to assess the performance of the entire method, we ran 
ANDROMEDA on real data from the VLT/NaCo instrument. We 
applied ANDROMEDA to data collected by the NaCo Large 
Programme, which aimed at dectecting planets on wide orbits 
(Chauvin et al. 2015). The data we used consisted in observing 
a bright star surrounded by numerous background stars that are 
natural point sources acting as companions to be detected by the 
algorithm. Synthetic companions were then added inside the im¬ 
ages in order to better quantify the ability of ANDROMEDA to 
accurately retrieve the position and contrast of the point sources 
present in the images. 

4.1. Data used for the test 

The data collected are sequences of saturated exposures (there 
is no coronagraph in that setup) taken in pupil tracking mode. 
The chosen star is TYC-8979-1683-1 (also called CD-62-657) 
observed in 2011 on May 5 within the ESO program 184.C- 
0567(D) (PI: J.-L. Beuzit). This star is a G7 star of 17 Myr 
(V = 9.36, H - 7.47) located at 75.6 pc from the Sun. The ob¬ 
servation was made in H-band (filter centered around 1. 66 / 1111 ) 
and stored in a 1024 x 1024 pixel frame (S13 camera inside the 
CONICA imager having a field of view of 13.6" x 13.6"), the 
pixel scale being 13.22 mas/pixel. 

4.1.1. Observation conditions 

The star was observed during a total integration time of 36 min 
(giving 11 cubes of 29 frames each with an exposure time of 
6.8 sec) and for a total field rotation of 18.5°. Seeing conditions 
were good but variable (seeing of 0.57"to 1.15"; coherence time 
of 4-9 ms; Strehl ratio of the reference PSPs: 21% and 26%). The 
empirical PSP core PWHM is estimated to be of 4.75 + 0.05 pix¬ 
els. The target was observed close to meridian crossing, the PSP 
core is saturated inside a radius of 10-15 pixels (0.13"-0.20"), 
and integration times are set short enough so that the angular 
smear of potential companions is small, especially in the inner 
part of the field. The data reduction of saturated exposures in¬ 
cluded sky subtraction (sky frame constructed from the median 
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Fig. 4. Images obtained by running the automatic detection module on the output provided by ANDROMEDA applied to the images of TYC- 
8979-1683-1 taken with the VLT/NaCo instrument (see Sect. 4). The threshold is set at 5cr. Left: S/N map subwindows (11x11 pixels) in which 
the 2D Gaussian fits are performed for each detection. The subimages are classified by decreasing S/N and indexed by an increasing integer from 
one to the total number of detections above threshold. If the fit could not converge, one asterisk is placed above the corresponding subimage. 
Otherwise, the two axes of the fitted Gaussian are superimposed and their lengths are equal to the FWHM. Right: S/N map showing the location 
of each detection in the first image (labeled by their index of detection) as well as their flux range (color of the index) and whether it is an artifact 
(one asterisk means that the 2D Gaussian fit did not converge; the @ symbol means that it is assessed as a tertiary lobe artifact). 


combination of exposures obtained at the five different jitter po¬ 
sitions: on minute timescales, the image center is moved by +3 
arcsec in x or y on the detector field. ), flat fielding, bad pixels 
correction, and rejection of poor-quality exposures. The location 
of the star center on each frame is determined by fitting the un¬ 
saturated portion of the saturated PSF with a 2D Moffat function. 
Individual frames are then shifted and registered to a common 
image center in between four pixels. 

A short series of exposures with the unsaturated star was 
taken before and after the main saturated sequence in order to 
build the reference PSF required as an input for ANDROMEDA. 
These unsaturated sequences were acquired with an exposure 
time of 1.7927 seconds in autojitter mode, using a neutral den¬ 
sity filter (ND_Short) of 1.19+0.05% transmission factor. The 
first unsaturated image obtained for these data is shown in Fig. 5- 
right. 

The parallactic angle associated with each frame of saturated 
images is computed from the observing time (converted from 
UTC to LST), assuming that individual exposures are recorded at 
constant time intervals within each data cube (time information 
is available only for the beginning and end of each data cube). 

We thus obtained the three necessary inputs for AN¬ 
DROMEDA: a reduced image cube, a PSE of reference, and 
a vector containing the parallactic angles of each image of the 
cube. 

4.1.2. Introduction of synthetic companions 

To better quantify and optimize the detection performance of 
ANDROMEDA using NaCo data, we implanted 20 additional 


synthetic substellar companions in the image cube. The signal 
of each synthetic companion was modeled using the unsaturated 
PSE image and inserted in the individual reduced frames, taking 
into account the field rotation that occurred between the expo¬ 
sures. The procedure used to inject synthetic planets in the im¬ 
ages of the data cube is explained in Chauvin et al. (2012). 

The twenty synthetic companions were introduced along five 
radial directions of respective position angles of 30°, 60°, 90°, 
120°, and 150° on the first image for each of these angle, at 
five angular separations of 0.26", 0.53", 1.06", and 2.12". The 
synthetic companions of the same position angle are of equal 
flux, each with peak intensities corresponding to magnitude dif¬ 
ferences of respectively 12.85, 12.10, 11.35, 10.60, and 9.85 for 
the five position angles in increasing order. An example of one 
saturated image of TYC-8979-1683-1 on which the location of 
these introduced synthetic companions are added is shown in 
Eig. 5-left. In order to verify the detection capability at close 
separation, we also added a companion at 0.26", with a PA of 
220° and a magnitude of 6.8 (contrast of about 2.10“^). 

4.2. ANDROMEDA’S output 

We applied ANDROMEDA to the data described above and ob¬ 
tained the four output detailed in Sect. 2. To process these data, 
we have chosen the user-defined parameters values gathered in 
Tab. 1 and whose suitability is discussed in detail in Sect. 4.2.3. 

The S/N map obtained is shown in Eig. 6 according to dif¬ 
ferent filtering fractions. It can be seen that when the low fre¬ 
quencies are removed, the companion candidates appear more 
distinctly, mostly those close to the star. The diffraction patterns 


Article number, page 9 of 19 






























A&A proofs: manuscript no. Andromeda 



-40 -22 -3 16 34 63 71 89 108 127 145(7 


Fig. 6. Effect of the pre-processing (spatial high-pass filtering of the raw data) on the output S/N map obtained by processing the TYC-8979-1683- 
1 NaCo data (including both real and injected synthetic companions) with ANDROMEDA. Left: S/N map obtained without pre-filtering. Middle: 
S/N map obtained with high-pass pre-filtering, using F = 1/16. Right: S/N map obtained with high-pass pre-filtering using F = 1/4. All these S/N 
maps are normalized following the procedure in Sect. 2.5 



Fig. 5. NaCo data of TYC-8979-1683-1. Top: One reduced saturated 
image of the raw data cube (600 x 600 pixels, linear scale). The position 
of the injected synthetic planets is indicated with colored dots on the 
green circles. Bottom: Reduced unsaturated image used as a reference 
PSE to perform ANDROMEDA (32 x 32 pixels). 


Table 1. User-defined parameter values chosen to run ANDROMEDA 
on the TYC-8979-1683-1 VLT/NaCo data. 


User-defined parameter 

Used value 

Filtering fraction, F 

1/4 

Minimum distance for subtraction, (Smin 

1.0 T/D 

Width of the annuli for ADI, dr 

1.0 T/D 

Ratio optimization to subtraction area, Ra 

2 

Size of the input PSF window, Npsj 

32 X 32 pixels 

Normalization profile smoothing, Nsmooth 

18 pixels 


not filtered enough, the number of false detections increases, es¬ 
pecially close to the star, and many companions that are close 
to the star and faint are missing. Moreover, the accuracy of the 
estimated position and flux decreases since the planetary pat¬ 
terns are wider and less regular. We consequently chose to set 
the filtering fraction to F = 0.25: a quarter of the low frequen¬ 
cies are removed. This value is a trade-off between not losing 
too much companion flux (see Sect. 2.1), detecting as many true 
companions as possible, and decreasing the overall number of 
false alarms. We note that the oversampling of this data set is of 
1.6, and so the filtering fraction must be smaller than F - 0.6. 

In Fig. 7, from left to right, are the flux map, the likelihood 
map, and the standard deviation of the flux map obtained with 
the parameters listed in Tab. 1. As expected, the flux and the 
likelihood maps look similar to the S/N map, and the standard 
deviation of the flux shows a clear radial decreasing trend from 
the center to the edge, as does the residual speckle noise distri¬ 
bution in the differential images. 


4.2.1. Effect of the normalization process 


induced by the spider of the telescope are revealed as two blurry 
lines crossing at the center (their cone shape is due to the ADI 
process per annulus and shows the expected pattern, negative- 
positive-negative); the patterns remain despite the ADI proce¬ 
dure, probably due to the smearing of the star in the field. If 


In order to correctly normalize the S/N map, we have plotted 
the radial profile of the robust standard deviation per annulus 
of the S/N map, according to different smoothing values. We 
then chose the best profile, the one that smoothed tiny irregu¬ 
larities and kept the global trend. Some profiles are presented 
in Fig. 8 to compare the robust versus regular standard devia- 
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Fig. 7. Output obtained by processing the TYC-8979-1683-1 NaCo data (including both real and injected synthetic companions) with AN¬ 
DROMEDA. Left: Flux map. Middle: Likelihood map. Right: Map of the standard deviation of the flux. For each of these maps (as for the 
S/N maps in Fig. 6) the central region corresponding to the star has been masked. The intensity scale is logarithmic for the likelihood and the 
standard deviation of the flux maps, and linear for the flux map. White corresponds to the highest value and black to the minimum value. 





Fig. 8. Azimuthal mean standard deviation of the S/N map as a function of the distance to the star TYC-8979-1683-1. Left: Regular radial profile 
of the S/N standard deviation. We can visualize the high peaks indicating the presence of companions at certain separations. Middle: Radial 
profile of the S/N robust standard deviation. The highest peaks have disappeared but the profile is still jagged. Right: Radial profile of the S/N 
robust standard deviation smoothed over 28 pixels to obtain the global trend of the S/N map standard deviation deprived of its companions. The 
ANDROMEDA process reduces the exploitable zone between the IWA (at 13 pixels here since IWA = 4A/D) and the OWA (at 281 pixels here 
since OWA = — Npgpll — dr). 


tion and smoothed versus unsmoothed profiles, thus justifying 
the choices that were made to build the normalized S/N map. 

The maps in Fig. 9 show the raw S/N map (left) and the 
normalized S/N map (right) both thresholded to 5cr. The non- 
normalized S/N map (left) does not provide any workable in¬ 
formation and just illustrates that its radial profile decreases to¬ 
ward the edge; however, the normalized S/N map reveals only 
the probable point sources in the map and therefore enables the 
automatic detection of companions. 

4.2.2. Interest of the filtering and normalization process 

In order to verify the consistency of the filtering and normaliza¬ 
tion procedure, we plotted the histograms of the residuals in the 
S/N map. This permits us to check if the PDF of the residuals 
has the expected shape, that is to say, Gaussian and white. 

The histograms shown in Fig. 10 have been produced from 
the normalized S/N map obtained by processing the filtered im¬ 
ages with ANDROMEDA on TYC-8979-1683-1, using the pa¬ 
rameters gathered in Tab. 1. It can be checked that the resid¬ 
ual noise in the S/N-map has been made Gaussian thanks to the 


high-pass filtering (Sect. 2.1) and also that, thanks to the normal¬ 
ization process, this Gaussian distribution is centered on zero, 
which means that the residuals have a mean value of zero, and 
its FWHM is equal to two, which corresponds to a standard de¬ 
viation of about one. 

4.2.3. Impact of the user-defined parameters 

This part is a short discussion of the influence of each of the 
user-defined parameters on ANDROMEDA’S output. The opti¬ 
mal values of the user-defined parameters found from these in¬ 
vestigations are set as default in the algorithm, though it is ap¬ 
propriate only for this particular TYC-8979-1683-1 data set. The 
parameters are discussed in the order they are used in the algo¬ 
rithm. We note here that most user-defined parameters do not 
have a major impact on the estimation performance. Their opti¬ 
mal values/or this data set, as well as the strength of their impact 
in terms of detection are gathered in Tab. 2. 

Performing ADI: To perform the ADI according to the tech¬ 
nique explained in Sect. 2.2.1, the parameter 6 mm must be cho- 
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Fig. 9. S/N maps that have been thresholded to 5cr. Left: Non-normalized thresholded S/N map. We observe that many signals are above threshold 
(white), mostly close to the center and along the spider diffraction pattern. Right: Normalized thresholded S/N map. Here only the probable 
planetary signals are found above the threshold. It is possible to perform an automatic detection on this map. 
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Fig. 10. Histograms of the residuals in the normalized S/N map, obtained using the NaCo data of TYC-8979-1683-1, inside four different annuli 
centered on the star: black solid lines are for an inner radius of 15 pixels (~ 5/1/D), purple solid lines for 50 pixels (~ 17/1/D), dark blue solid lines 
for 95 pixels (~ 32/1/D), and light blue solid lines for 247 pixels (~ 90/1/D). Each annulus has a width of 15 pixels except for the largest, which 
only has a width of 8 pixels. No obvious planetary-like signal can be found inside these annuli. Left: without filtering. Middle: with a pre-filtering, 
using F = 1/16. Right: with a pre-filtering, using F = 1/4. Gaussian fits of the histograms are overplotted in red dashed lines. 


sen carefully to avoid self-subtraction of the companion while 
subtracting images as temporally correlated as possible. For this 
data set, the best compromise is obtained when using = 
l.O/l/D. 

Since the displacement between two given images varies 
with the separation from the star, the ADI is performed over 
an annulus of user-defined width dr that should be the thinnest. 
Indeed, there are slightly fewer artifacts detected when dr de¬ 
creases. However, the width of the annuli should not be too small 
because it increases the processing time. A compromise is to set 
dr — l/l/D as default. 

When choosing to perform the image subtraction with ei¬ 
ther a least-squares fit or a LI-affine fit, the width of the annulus 


over which is calculated has a certain ratio Ra with respect to 
the width dr of the subtracted annuli (see Sect 2.2.1). There is 
a compromise to make when choosing the optimization area so 
that it can avoid discontinuities (being thicker) while minimizing 
the flux difference within the annulus taken into account (being 
close to the subtraction area). As expected, when is too large 
many artifacts may be detected, but when Ra is too close to one 
there are discontinuities appearing between adjacent annuli. Un¬ 
der these considerations, the optimal value is set to Ra - 2. 

Maximum likelihood: The size of the PSF window, Npsj, must 
completely enclose the full signal, otherwise it induces strong 
errors on the flux estimation. On the other hand, the larger the 
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window, the longer it takes to run ANDROMEDA. Thus, this 
parameter can be optimally set so that the first bright ring fits 
inside the PSF window. 

We note that the four parameters talked about previously 
i^min, dr, Ra, and Npsf) are parameters that have been pointed 
out from the first developments of the ADI method. Some algo¬ 
rithms such as TLOCI (Marois et al. 2014) or PC A (Amara & 
Quanz 2012; Soummer et al. 2012), which are improvements of 
ADI, try to reduce the number of parameters for this part of the 
algorithm. 

Normalization of the S/N: The number of pixels on which the 
S/N robust radial standard deviation profile is smoothed, Nsmooth, 
is another important parameter. Too much smoothing biases the 
normalization and provokes the appearance of artifacts detected 
at less than 0.5" from the star but that can be easily rejected 
a posteriori. On the other hand, a low smoothing may result in 
missing faint close signals. Thus, a good value is to smooth the 
profile over about 28 pixels, but this value is strongly dependent 
upon the oversampling and upon the quality of the images. It is 
the only parameter that must currently be chosen by hand after 
visualizing the trend of the profile with different values. We note 
that the chosen value does not need to be very precise: for this 
data set, it can vary from about 26 to about 50 pixels without 
adding too many artifacts above threshold. 

To conclude, along with a correct normalization, AN¬ 
DROMEDA provides workable output that allow an automatic 
detection procedure to be built. The results obtained with the au¬ 
tomatic procedure, developed to extract the position and flux of 
the candidate companions from the output maps, are presented 
in the next section. 

4.3. Analysis of the planetary signals detected 

In a second step, the automatic detection module was applied 
on the maps, using a threshold of 5cr. The subimages extracted 
from the S/N map in which the fits are performed are shown in 
Fig. 4-left. On account of the pixel scale of the imaging camera, 
the size of the subwindows is set to 11 pixels, which is the ex¬ 
pected planetary pattern size (~ 4/1/D) that encloses the whole 
planetary signal. The S/N map on which the position of the de¬ 
tected signals is visible and labeled by their index of detection is 
shown in Fig. 4-right. 

The program detected 39 signals above this threshold, in¬ 
cluding 25 reliable detections (Sect. 3.3). Of course the number 
of detections found in the S/N map depends on the threshold 
set; this dependence is dealt with in detail in Sect. 4.4. As ex¬ 
pected, only the highest S/N signals are found surrounded by 
tertiary lobes that are detected (detections #1, #2, and #3, having 
S/N values of respectively 150(T, 105cr, and 95cr). It is also no¬ 
ticeable that, as expected again, all the fits that have converged 
in the S/N map have also converged in the flux map and vice 
versa. Thanks to the criteria classiying the detections laid down 
in Sect. 3.3 , the automatic detection procedure efficiently sepa¬ 
rates probable true companions from artifacts. 

4.3.1. Performance in terms of detection and analysis 

To better quantify for the performance of the method in terms of 
detecting point sources present in the image field and in estimat¬ 
ing their positions and flux, it is possible to use the knowledge 
we have about the injected synthetic planets. Knowing their ex¬ 
act position and contrast, it is possible to compare them with the 


values obtained by the algorithm. It is also important to check 
that all the planetary companions found are indeed above the de¬ 
tection limit derived from this data set. 

Figure 11 shows the contrast in terms of magnitude of the 
detected point sources as a function of the angular separation be¬ 
tween the detected signal and the central star. On this graph, six 
horizontal lines represent the original contrast of the synthetic 
companions and four vertical lines are placed at their theoreti¬ 
cal separation from the star. Consequently, we know that signals 
from the synthetic planets ought to be found at every crossing 
between the horizontal and vertical lines, except for the upper 
horizontal line, which only has one synthetic companion at the 
closest separation. The detection limit is overplotted on the graph 
(solid line) and the error bars in terms of 3(T error on the contrast 
estimation are added to the graph. 

We first notice that, apart from the brightest synthetic com¬ 
panion that we injected just above the detection limit, none of 
the five synthetic planets that had been injected at 0.26" from 
the star is detected. Moreover, only the brightest synthetic planet 
located at 0.53"is detected. This result is expected since each of 
these undetected synthetic planets is located under the detection 
limit curve. The detection limit shows the following trend: close 
to the star only bright companions can be detected, whereas go¬ 
ing farther from the star enables fainter signals to be detected. Of 
course, the detection limit is dependent upon the chosen thresh¬ 
old, but even with a very low threshold, it is impossible to detect 
these faint and close signals. As all image processing methods, 
ANDROMEDA is limited by the observation conditions. 

The estimated positions of the detected synthetic planets 
match the theoretical values, even very close to the star. The er¬ 
ror bars in position due to the Gaussian fit uncertainty at 3cr are 
of about 2.0 mas (from 0.02 mas for the highest signal to 6.0 
mas for the faintest). These errors in position are not shown in 
Fig. 11 since they are smaller than the symbol size. As expected, 
the flux is better estimated when the companion is bright and far 
from the star. A good agreement is still observed between the 
theoretical contrast and the estimated value, knowing that the er¬ 
ror bars shown in Fig. 11 are only the ones given by the map of 
the flux standard deviation but without taking into account the 
instrumental errors or the algorithm’s intrinsic errors. 

To assess the errors intrinsic to the algorithm, one approach 
is to slightly move the user-defined parameters (which might in¬ 
fluence the estimations) from their optimal value. The bound¬ 
aries within which such user-defined parameter are made to 
vary are the following: d,,,,-,, 6 [0.2; 2.0]/l/D, dr e [0.5;3]/i/D, 
Ra e [1;4], and Nsmooth must be varied experimentally from no 
smoothing to what can be reached with the data at hand. In this 
way, we obtained an error of about 5.0 mas in position and of 
about 0.25 magnitude, depending on the intensity of the signal. 
Accounting for this deviation, the real position and contrast val¬ 
ues of the injected fake companions are within the error bars of 
the ANDROMEDA estimations. 

In exoplanet imaging, the dominant errors on the position 
and flux estimations are usually due to instrumental sources. But 
in the case of poor quality data, it must be verified which of the 
instrument calibration errors or the algorithm unstability errors 
are dominant. 

4.4. Threshold sensitivity 

This section discusses the optimal threshold range that would 
reveal as many true companions as possible, while not miss¬ 
ing any. In particular, it is tested whether a constant detection 
threshold is efficient all over the field of view. This discussion re- 
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Table 2. User-defined parameters set as defaults in the ANDROMEDA pipeline and their respective significance. 


Parameter 

Definition 

Units 

Default value 

Impact 

F 

Filtering fraction (Sect. 2.1) 

- 

1/4 

low 

^min 

Minimum separation to build the differential images (Sect. 2.2.1) 

A[D 

1.0 

high 

dr 

Width of annuli on which ADI is performed (Sect. 2.2.1) 

A[D 

1 

low 

Ra 

Ratio optimization to subtraction areas (Sect. 2.2.1) 

- 

2 

low 

Npsf 

Size of the square PSF image (H-band filter of NaCo) 

pixels 

32 

see text 

^smooth 

Smoothing of the S/N robust standard deviation profile (Sect. 2.5) 

pixels 

18 

high 


Notes. The right column shows the impact of the user-defined parameters over the whole process: If low, the value can remain fixed and if high, it 
should be tuned according to the data set. 
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Fig. 11. Contrast of the detections as a function of their distance to the star TYC-8979-1683-1. The theoretical contrasts and distances of the 20 
injected synthetic planets are shown as vertical and horizontal lines: we should have found a planetary signal at every crossing. One synthetic 
companion has also been added at the closest separation just above the detection limit in order to check its consistency (red dashed line). The 
detection limit at Scr is overplotted (solid line). This detection limit curve is corrected for the underestimation of the small sample statistics at 
close separation from the star (see Sect. 2.4). The minimum value of the flux standard deviation map at 5(T for each separation is overplotted for 
information (dotted line). The radial median of the flux standard deviation map at 5(T, however, provides a realistic estimation of the detection 
limit. The detected signals assessed to be tertiary lobes are not shown on the graph. 


lies upon a study of the behavior of the number of detections as a 
function of the threshold. Discussing missed detections and false 
positives requires assumptions on the actual number of compan¬ 
ions present in the field. The series of fake companions is a firm 
basis. The number of astronomical background stars produces 
additional test cases to check the detection capability homogene¬ 
ity, and can be tested in comparison to other ADI approaches. 
We assume in the following that the exact number of detectable 
point sources for this data set is 25. The graph in Fig. 12 was 
obtained by running ANDROMEDA under the same conditions 
as before (Tab. 1), and represents the number of detections as a 
function of the threshold which varies from 3cr to 6cr. On this 
graph, both the total number of detections and the number of 
so-called reliable detections are indicated. 

This graph shows the existence of a short optimal threshold 
range, between 4.7cr and S.lcr, for which the number of detec¬ 


tions is exactly the one suspected of being the true one. When in¬ 
creasing the threshold above this range, the faintest signals (like 
#39) and/or the ones very close to the star (like #33) are the first 
to be missed by the automatic detection process, whereas sev¬ 
eral tests have proved that these are most likely to be real point 
source signals. When decreasing the threshold, some detections 
that are actually not true signals are detected and still pass the 
test that is supposed to sort out artifacts from potential real sig¬ 
nals. These detections are usually faint signals (AH ~ 14-16), 
located far from the star where the algorithm should show bet¬ 
ter performance, that are most likely to be false alarms due to 
residual speckles. Moreover, a careful examination of the shape 
of this graph brings useful insights to help the user choose the 
threshold. An initial inspection shows that lowering the thresh¬ 
old increases the number of detections exponentially and that 
there is a plateau, from ~ 4.5cr to ~ 5.0cr, where the number 
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Fig. 12. Number of total detections (dark color) and of confirmed detec¬ 
tions (light color) as a function of the threshold set for the TYC-8979- 
1683-1 NaCo images. The optimal zone (where all 25 companions are 
detected) is shown in green. A false alarm zone is shown in red (more 
than one false detection) and orange (one false detection) and a loss 
zone in blue (more than one true companion is not detected) and purple 
(one true companion is not detected). 

of detections is quite stable as it is for the number of true de¬ 
tections. Also, by looking at the ratio of true detections versus 
ill-fitted ones, below a threshold of 4.5cr, the number of ill-fitted 
detections becomes significantly higher than the number of true 
detections, whereas above the threshold this ratio remains quite 
constant. This study is of course strictly valid only for this data 
set, which has the advantage of containing many point sources, 
but it should still give a fair idea of the behavior of the algorithm 
as a function of the threshold. 

Even if this test case is very advantageous, the goal of this 
section is to check the overall method skills. This procedure 
proved to be very efficient in terms of detection ability and in 
terms of position and contrast estimations with low error bars. 
We noticed that thresholding the S/N map was not enough to 
provide only real probable planetary signals so we added an 
a posteriori classification of the detections, which also proved 
its efficiency in discriminating artifacts from very probable true 
point source signals. We note that ANDROMEDA does not pro¬ 
vide more information than can be reached given the observa¬ 
tion conditions (according to the detection limit trend), but the 
algorithm does provide, within some minutes, reliable informa¬ 
tion in terms of position and flux of the detected companions. 
Of course where the noise is higher (toward the star), the esti¬ 
mation is slightly less accurate, but the errors still remain under 
the observational errors such as the PSF centering or the tip-tilt 
correction. Section 5 is a comparison of the results obtained by 
processing the well-known case of Beta-Pictoris using the PCA 
methods or ANDROMEDA. 


5. Results on p Pictoris and comparison with the 
PCA-KLIP method 

In addition to testing the method on the previous field, popu¬ 
lated with both numerous background stars and additional syn¬ 
thetic companions, we also applied it to the emblematic and 
well-known case of (3 Pictoris. This star is surrounded by a de¬ 


bris disk inside which only one planet, p Pictoris b, has been 
detected by imaging. This close companion is located at ~ 9 
AU from its host star and was first detected by Lagrange et al. 
(2009). This companion has since been observed repeatedly be¬ 
cause of its outstanding astronomical interest for many reasons. 
It is a moderate-mass giant planet at close physical separation, 
making it a good candidate for having been formed within the 
disk. Its interactions with the remaining debris disk can still be 
witnessed and dynamical constraints can also be derived from 
the orbital motion and from numerous falling evoporating bod¬ 
ies as observed in spectroscopy (see Lagrange-Henri et al. 1988; 
Beust et al. 1990). In the context of this paper, observations of 
this companion provide an excellent test case of our method in 
the challenging case of the detection of a companion at very 
short separation (< 0.5") with appropriate high-quality corona- 
graphic images; the images processed by ANDROMEDA should 
be compared to those processed thanks to other methods repre¬ 
senting the state of the art. 

The p Pictoris data on which we applied ANDROMEDA are 
described in Absil et al. (2013). The goal of this data set was to 
search for closer planetary companions (down to 2 AU from the 
star), which would be at the origin of the dynamical excitation in 
the planetary system. This search required the use of a newly de¬ 
veloped coronagraph that allows detections at angular separation 
as close as O'.T. 


5.1. Data used for comparison 

P Pictoris was observed on 31 January 2013 for 3.5h at the Very 
Large Telescope (Chile) using the NaCo instrument (ESO pro¬ 
gramme ID 60.A-9800). The images were recorded in the L' 
band (centered on 3.8jum), which is commonly used to work 
in a more favorable planet/star contrast regime compared to 
shorter wavelengths, and to obtain an image quality close to 
that delivered by extreme adaptive optics (XAO) systems in the 
near-infared. To minimize the impact of starlight in the images, 
the newly commissioned L'-band Annular Groove Phase Mask 
(AGPM, Mawet et al. 2005, 2013) vector vortex coronagraph 
was used. This kind of coronagraph provides an inner work¬ 
ing angle and an achromaticity over the bandwidth of interest 
as good as any other phase mask developed so far. 

The data are constituted of ten blocks of 200 successive 
frames of 0.2s, each taken under fair turbulence conditions. The 
observing sequence was obtained in pupil-tracking mode, show¬ 
ing a total field rotation of 83°. Two non-coronagraphic PSFs 
were taken before and after the observing sequence by moving 
the star far from the vortex center. The pixel scale of the camera 
used is 27.15 mas/px and the PSF FWHM is empirically mea¬ 
sured to be 4 pixels. All the information concerning the data set 
used in this section and the pre-processing steps applied before 
running image processing algorithms such as ANDROMEDA 
(basic cosmetic treatment, recentering, frame selection, etc.) can 
be found in Absil et al. (2013). 

5.2. ANDROMEDA results and comparison with PCA 

This section compares the quality of the detection, photometry, 
and astrometry obtained by ANDROMEDA with respect to the 
output of a principal component analysis (PCA Soummer et al. 
2012; Amara & Quanz 2012), as presented in Absil et al. (2013). 
We thus ran ANDROMEDA on the exact same data cube and the 
automatic detection module on its output. 
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Fig. 13. Resulting S/N map of jj Pictoris according to two different image processing methods applied on the reduced data. Left: S/N map obtained 
from ANDROMEDA, using a least-squares optimization for the ADI subtraction. Middle: S/N map obtained from ANDROMEDA, using a Ll- 
affine optimization for the ADI subtraction. The residuals are slightly lower, mostly close to the central part, as explained in Sect. 2.2.1. Right: 
S/N map obtained with the PCA-KLIP algorithm (see text). All images are linearly scaled. No signal is found above the 5cr threshold in any of the 
maps, except for the companion f Pictoris b. 

Table 3. Relative astrometry and photometry of f Pictoris b retrieved from the VLT/NaCo-AGPM data processed with PCA-KLIP and AN¬ 
DROMEDA methods. 


Parameters 

Estimated value using ANDROMEDA 

Estimated value using PCA 

Peak S/N (cr): 

17 

17 

Separation from the star sep (mas): 

455 + 8.7 

452 + 9.6 

Position angle PA (°): 

210.7 + 0.8 

211.2+1.3 

Contrast AL' (mag): 

8.09 + 0.21 

8.01+0.16 


Given that for this study we have selected a region of 
150 X 150 pixels around the star, and knowing the pixel scale 
of the L' camera (27.15 mas/pixel), the user-defined parameter 
values chosen to run ANDROMEDA and the detection module 
are listed in Tab. 4. 

Table 4. Parameters used to run ANDROMEDA on the f Pictoris VLT- 
NaCo-AGPM data set. 


Parameter 

Used value 

F 

1/4 (default) 

dr 

1 T/D (default) 

Ra 

2 (default) 

^min 

1 T/D (default) 

Npsf 

20 X 20 pixels 

^smooth 

10 pixels 

IWA 

1 T/D 

Threshold 

5cr 

Subwindow size 

9x9 pixels 


The resulting S/N map is shown in Fig. 13-Left where the ex¬ 
pected companion f Pictoris b is clearly visible as a sharp bright 
spot, southeast of the star (located at the center of the frame). 

The same data set has been processed with a home-grown 
PCA algorithm, based on the KLIP implementation (Soummer 
et al. 2012). The KLIP algorithm uses a truncated basis of eigen¬ 
vectors created by a Karhunen-Loeve transform of the initial set 
of images, to perform the subtraction of the star residuals. To fol¬ 
low the original KLIP algorithm, the image processing was per¬ 
formed on the full 150 x 150 pixel frames at once, and no frame 
was excluded from the ADI image library based on the amount 


of parallactic angle variation (unlike in Absil et al. 2013). From 
the final image of the PCA processing, we compute a S/N map 
by testing resolution elements centered on each pixel of the map 
against all the other resolution elements located at the same an¬ 
gular distance from the star, as described in Mawet et al. (2014). 
The number of principal components used in KLIP was tuned to 
maximize the S/N of the planet in the final image, resulting in 
the use of 18 principal components instead of 30 in Absil et al. 
(2013). The S/N map is illustrated in Fig. 13-Right, showing a 
peak S/N ~ 17 on the planet, i.e., the same as in the case of AN¬ 
DROMEDA. We note that more advanced implementations of 
the KLIP algorithm, working on well-defined subregions (e.g., 
annuli or parts of annuli) in the images and including a frame re¬ 
jection criterion based on the parallactic angle, can reach a peak 
S/N of up to ~ 20, although at the expense of a drastically in¬ 
creased computation time (a few minutes instead of ~ 1 sec¬ 
ond for the classical KLIP algorithm). Using ANDROMEDA, 
we could reach this S/N value by smoothing the profile over 20 
pixels instead of ten, as used here. This does not affect the com¬ 
putation time, but for this value of Nsmooth one artifact appears 
above the 5cr threshold. 

The errors on the photometric and astrometric estimations 
are mostly limited by instrumental calibration errors, which can 
be decomposed into three main contributions: the error on the 
position of the source (8.5 mas), the PSF centering error in the 
images (0.1 mas), and the plate scale error (0.04 mas). As these 
errors are independent, the resulting error is the quadratic sum, 
giving in this case a total of 8.5 mas. Adding the 3cr error on 
the position estimation with ANDROMEDA (1.7 mas), the total 
error on the position is of 8.6 mas and the algorithm’s intrinsic 
error is evaluated within this error bar. For the position angle, 
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the main errors are due to the true north direction knowledge 
(0.09°) and the offset of the derotator (0.01°), as discussed by 
Chauvin et al. (2012). ANDROMEDA gives the PA with a 3cr 
precision of 0.8° and a negligible intrinsic error. The errors on 
the flux estimation are mainly due to the variation of the PSF 
along the observation (0.05 mag). ANDROMEDA estimates the 
flux with a 3cr statistical error of 0.2 mag and a negligible in¬ 
trinsic error, and as these three sources are independent again, 
the total error is the quadratic sum of each error. The estimated 
position and contrast with both methods and their corresponding 
error bars are found in Tab. 3, where the photometric and astro- 
metric estimations for the PCA pipeline are directly taken from 
the publication of Absil et al. (2013). The two estimations agree 
with each other within error bars. 

We emphasize once again here that the photometry and 
astrometry can be recovered much more easily with AN¬ 
DROMEDA than with PCA, in particular, because the estimation 
is direct and precise without needing to inject fake companions. 
The same applies for the detection limit computation, which is 
given directly by the map of the standard deviation of the flux 
provided by ANDROMEDA and does not need, for instance, nu¬ 
merous synthetic planet injections as for other methods. 

5.3. Sensitivity comparison and discussion 

Another way of comparing the performance of the two methods 
is to plot their 5cr detectability limits in terms of contrast. 

Eor the PCA-KLIP method, the computation of the detec¬ 
tion limits is based on the hypothesis that the noise is Gaussian 
(Marois et al. 2008), which was checked empirically in Absil 
et al. (2013). Several methods have been used in the literature 
to derive detection limits. Here, we use the statistical framework 
presented in Mawet et al. (2014), where test speckles are com¬ 
pared to the statistics of the flux within all other resolution ele¬ 
ments located at the same angular distance from the center. This 
framework assumes that the noise statistics varies with the dis¬ 
tance from the star but not with the azimuth, and it takes into 
account the effect of small sample statistics on the noise estima¬ 
tions. More specifically, the computation of the contrast curve 
consists in normalizing the flux in the final PCA-processed im¬ 
age by the integrated flux contained inside the (off-axis) stel¬ 
lar PSF core, then in calculating the standard deviation of the 
fluxes enclosed in all the non-overlapping apertures of diameter 
equal to the FWHM of the PSF that can be placed at a given 
distance from the center. The same procedure is repeated at all 
radial distances. During the process, any known companion can 
be masked out, which we did here in order not to bias our noise 
estimation at the distance of fi Pic b. This estimation of the con¬ 
trast curve provides more realistic results than the pixel-to-pixel 
standard deviation noise estimation that has frequently been used 
in the literature in the past. 

We chose to use here a “smart” version of the PCA-KLIP al¬ 
gorithm (sPCA), where the images are decomposed in annuli and 
where a frame rejection criterion is included based on the par¬ 
allactic angle, in order to avoid self-subtraction of the possible 
planetary signals. We used the same smart PCA parameters as 
for the computation of the detection limits in Absil et al. (2013). 
The final contrast curve, plotted in Fig. 14, takes into account 
the algorithm throughput, estimated at several radial distances 
(and azimuths) by injecting fake companions in the original data 
cube, running the PCA pipeline with the exact same parameters, 
and comparing the flux of the fake companions in the final PCA- 
processed image with the flux of the injected companions. We 
note that the PCA contrast curve presented here is a factor of 


two higher than the one presented in Absil et al. (2013), where 
an improper convolution by a Gaussian kernel led to a factor of 
two underestimation of the contrast curve. 

Figure 14 shows the different detection limit curves obtained 
with the sPCA algorithm and with ANDROMEDA using the 
methods noted above, all corrected for the small sample statis¬ 
tics. The location of p Pictoris b in terms of angular separation 
and contrast is also indicated on the graph for comparison. In 
this figure, we can see that at angular separations ranging from 
0.1 "to 2", ANDROMEDA and PCA provide similar detection 
limits. 

To give a rough idea, we mention here that ANDROMEDA 
takes about 180 seconds to process the 612 150 x 150 pixels im¬ 
ages on a quadri-pro six-core Intel Xeon. This computing time 
includes the time to obtain the three main outputs from AN¬ 
DROMEDA (estimated flux map, S/N map, and the map of the 
standard deviation of the estimated flux) plus the detection de¬ 
scribed in Sect. 3 (detection of probable point source, sorting out, 
subpixel position estimation, flux estimation, and corresponding 
derivation of the 3cr error bars and the computation of the de¬ 
tection limit curve). This latter detection is obtained from the 
ANDROMEDA output within a few seconds. For the PCA-KLIP 
method, it takes about 1 second to produce the final reduced im¬ 
age using the original KLIP algorithm. Building the S/N map in 
Fig. 13 from this reduced image along the Mawet et al. (2014) 
prescription takes an additional few tens of seconds. We note 
that the annulus-wise, “smart” PCA (sPCA) processing used to 
compute the detection limits in Sect. 5.3 takes about 10 min¬ 
utes on a recent core-i7 laptop. Also, we emphasize here that 
the photometry and astrometry are included in ANDROMEDA 
while PCA needs dedicated, computationally intensive methods 
to get an accurate estimation (such as the negative fake compan¬ 
ion technique, which requires running many PCAs to obtain the 
photometric and astrometric estimations together with their er¬ 
ror bars). The latter figures must be taken with care since there 
are still ways to optimize the computation time for both algo¬ 
rithms (for instance parallelizing the ADI plus MLE with AN¬ 
DROMEDA), but they are mentioned here to give a rough order 
of magnitude of the processing time. 

The two algorithms that we have used to process the 
VLT/NaCo-AGPM data of yS Pictoris both efficiently detect the 
close companion yS Pictoris b. However, ANDROMEDA does 
the detection automatically, by giving direct access to the S/N 
map (which is only a by-product in the PCA algorithm). The ap¬ 
proaches that the two algorithms use to estimate the position and 
flux of the companion are quite different, but result in compati¬ 
ble values. ANDROMEDA potentially shows better accuracy in 
terms of separation and contrast estimation since it directly re¬ 
lies on the detection probability derived from a strong hypothesis 
about the image formation and the residual noise in the images, 
which are well verified. To conclude, the two methods are com¬ 
plementary since they do not rely on the same inner concept and 
it is better to use several methods that are significantly different 
to better judge of the truthfulness of a detection. 

6. Conclusion and perspectives 

The ANDROMEDA method is designed to detect planetary-like 
signals in high-resolution and high-contrast images. In this pa¬ 
per we have described the improvements brought to the orig¬ 
inal method when applied to on-sky VLT/NaCo data, in order 
to make it more robust. This included the implementation of a 
pre-processing (high-pass filtering of the data), a post-processing 
(normalization of the S/N map), and the use of a robust method 
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Angular separation from star (mas) 

Fig. 14. Detection limits at 5cr obtained by processing the p Pictoris data from VLT/NaCO-AGPM with ANDROMEDA (green solid line, when 
using a Ll-affine optimization or blue dashed line when using a LS optimization) or sPCA (red dash-dotted line). The asterisk symbol shows the 
location of the detected signal above the 5<r threshold set (only p Pictoris b is detected here). 


to perform the image difference. As a result, this method proves 
to be highly efficient at processing experimental data from the 
VLT/NaCo instrument taken in pupil-tracking mode. In partic¬ 
ular the produced normalized S/N map can be thresholded with 
a constant value throughout the field of view so as to provide 
detectability and false alarm levels that are homogeneous over 
the field. This allows one to automatically detect the candidate 
companions present in the image, and classify them according 
to their S/N values. Artifacts can then be rejected based on mor- 
phogical criteria expected for a true companion. We then showed 
on synthetic companions and on the well-known case of p Pic b 
that ANDROMEDA accurately estimates the subpixel position 
and the contrast of these companions. Notably, and unlike most 
other methods, the selection of the data reduction parameters 
(like filtering or minimum angular separation used for image dif¬ 
ferences) is fully taken into account in the search for the com¬ 
panion signature, and thus it does not bias the companion flux 
estimate. In particular, the analysis of the result does not require 
subsequent tests with synthetic planet injections to estimate a 
companion flux loss. 

The second generation of instruments dedicated to exoplan¬ 
etary detection and characterization, such as VLT-SPHERE, 
(Beuzit et al. 2008), Gemini-GPI, (Macintosh et al. 2008), and 
Subaru-SCeXAO, (Martinache & Guyon 2009; Jovanovic et al. 
2013), is about to deliver a large amount of data that will re¬ 
quire massive, homogenous, and efficient companion extraction. 
The capabilities of ANDROMEDA in terms of performance and 
automatic detection will be very beneficial in this context. 

Though ANDROMEDA is operational in its current state, 
any further experience on these new instruments may motivate 
the evolution of this algorithm^ In particular, its probabilistic 
formalism directly enables one to include any a priori knowl- 


* The authors of the ANDROMEDA code are open to working in col¬ 
laboration in order to apply the code to other data and thus gain further 
experience. 


edge of the noise structure of the data or on the sought com¬ 
panions. Evolutions can be considered in various ways to build 
differential images (for instance, if it appears that the closest im¬ 
ages in time are not necessarily the most correlated). Einally, the 
non-coronagraphic PSE, defining the companion’s expected sig¬ 
nature, is currently assumed to be stable along the data sequence, 
as measured before or after the coronagraphic or saturated data 
set. If the evolution of the PSE and in particular the evolution of 
the Strehl ratio are known, this variability can be easily included 
in the algorithm for a better flux estimate fidelity. 
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